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Abstract 

A new family of parameters intended for composition studies in cosmic ray sur- 
face array detectors is proposed. The application of this technique to different array 
layout designs has been analyzed. The parameters make exclusive use of surface 
data combining the information from the total signal at each triggered detector and 
the array geometry. They are sensitive to the combined effects of the different muon 
and electromagnetic components on the lateral distribution function of proton and 
iron initiated showers at any given primary energy. Analytical and numerical stud- 
ies have been performed in order to assess the reliability, stability and optimization 
of these parameters. Experimental uncertainties, the underestimation of the muon 
component in the shower simulation codes, intrinsic fluctuations and reconstruction 
errors are considered and discussed in a quantitative way. The potential discrimina- 
tion power of these parameters, under realistic experimental conditions, is compared 
on a simplified, albeit quantitative way, with that expected from other surface and 
fluorescence estimators. 
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1 Introduction 



Ultra-high energy cosmic rays (UHECR) are studied by the detection of the 
extensive air showers (EAS) that they produce in the atmosphere. At present, 
there are two main observation techniques. One, the so called fluorescence 
technique, employs telescopes to collect the fluorescence light emitted by at- 
mospheric Nitrogen molecules after being excited by the charged particles of 
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the cascade. The other one, the surface technique, is based on the analysis 
of a discrete sampling of the shower front at the ground level, performed by 
surface detectors. UHECRs produce few observables. They are, basically, the 
arrival direction, the energy and some statistical hint about the identity of 
the primary particle. Of these three pieces of information, the geometrical one 
is the most reliable for both detection techniques [1,2,3,4,5]. The energy of 
the shower can be inferred with an accuracy of around 20% in the case of 
stereo fluorescence reconstruction [6] . In the case of hybrid detectors, this ac- 
curacy can be transferred to surface arrays using as calibrators events which 
are observed simultaneously by both techniques. The Auger Observatory was 
the pioneer of this strategy [7] which is likely to become a standard for accu- 
racy in the field [5,8]. The need for such cross calibration already highlights 
the existence of as yet either unidentified problems with our understanding of 
the physics involved in shower generation and development or inconsistencies 
in our implementation of those physical processes into the available shower 
simulation codes. Most of these problems certainly have their roots in the 
uncertainties associated with the extrapolation of cross sections, multiplici- 
ties and inelasticities from accelerator measurements at much lower energies, 
required to treat the first hadronic interactions experienced by the incoming 
cosmic ray in the upper layers of the atmosphere. Given the indirect nature 
of the detection of cosmic rays at the highest energies, those uncertainties 
permeate, to a larger or lesser extent, all measurements done afterwords. In 
particular, they have their strongest manifestation on the interpretation and 
reliability of mass composition tracers, since variations in cross section or in- 
elasticity can easily be misinterpreted as changes in baryonic composition (see 
[9] and [10] for a thorough discussion on the relevance of hadronic interaction 
models to shower development). 

Each observation technique has specific composition indicators. The fluores- 
cence technique is, at present, the most reliable one for composition studies 
since the longitudinal development of the charged component of the atmo- 
spheric shower is a relatively straightforward observable. Differences in com- 
position, manifest themselves through differences in the cross section for in- 
teractions with atmospheric nuclei. These, in turn, are mapped as different 
average depths of maximum development of the electromagnetic component 
in the atmosphere {X max ) for different nuclei and as different statistical disper- 
sions for the X max distributions of different nuclei (AX max ). If, for example, 
proton and Iron primaries are compared, the smaller cross section of the for- 
mer will produce larger X max and AX max than the latter [11,12]. However, 
unforeseen changes in cross section as a function of energy can affect these 
parameters in much the same way as true changes in composition would [13]. 

Surface detectors, on the other hand, sample the lateral distribution function 
(LDF) of the EAS at discrete points while they traverse the ground level. 
Beyond a few tens of meters from the shower axis, the particle content of the 
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shower at ground level is dominated by just two components, electromagnetic 
(i.e., electrons, positrons and photons) and muonic. These two sets of particles 
propagate differently through the atmosphere: the electromagnetic component 
propagates diffusively, while the muons do so radially from the last hadronic 
interaction region, which produces their parent mesons. Therefore, in a simpli- 
fied way, the shower front can be thought of as the composition of two fronts, 
a muonic one, which arrives first and is temporally thin and an electromag- 
netic disk, more extended in time that follows the muon front. Furthermore, 
the muon shower front has a much better defined and larger curvature radius 
than the electromagnetic front. One of the practical effects of these differences 
is that information about the relative intensity of both components inside the 
shower and, therefore, of the identity of the primary nuclei, is encoded simulta- 
neously in a non-trivial way in several properties of the shower front at ground 
level: the slope of the lateral distribution function of particles, the radius of 
curvature and its time structure. Thus, several parameters have been proposed 
to extract composition information from the surface measurements of EAS. 
The discrimination power of each parameter depends on the type of detector 
selected, water Cherenkov tanks (e.g., Haverah Park, Auger) or scintillators 
(e.g., Volcano Ranch, AGASA, Telescope Array). The former are sensitive to 
both the electromagnetic and muonic components (although the number of 
muons is much smaller than the electromagnetic particles, they generate a 
large amount of Cherenkov photons in the tanks) while the latter are sensitive 
to the charged particles of the cascade which are dominated by electrons and 
positrons. Scintillators are also shielded or buried underground to measured 
the muonic component of the showers. The most useful parameters so far have 
been the slope of the LDF [14,15], the curvature of the shower front, the num- 
ber of muons in the shower [16,17], several indicators of the time structure of 
the shower like the rise time and other related parameters [14,18,19], and the 
azimuthal asymmetries in the signals [20]. 

In general terms, fluorescence composition indicators are regarded as easier to 
observe and interpret, as well as less prone to systematic errors than surface 
parameters do. However, fluorescence detectors are constrained by a duty cycle 
of approximately 10%, while well operated, stable surface detectors can reach 
duty cycles near 100%. This factor alone, which makes the statistics per unit 
time of surface arrays an order of magnitude larger than that of fluorescence 
detectors, gives a great attractive to the search for reliable surface composition 
parameters. 

In the present work, we propose a new surface parameter which, we argue, for 
the same integration time can deliver better discrimination power than X max . 
The proposed parameter is defined as: 




(1) 
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where the sum extends over all the triggered stations N, tq = 1000 m is a 
reference distance, Si is the signal measured at the ith station and rj is the 
distance of this station to the shower axis in meters. Sf, is sensitive to the 
combined effects of the different muon and electromagnetic components on 
the lateral distribution function of proton and Iron initiated showers. 

As a case study, we analyze hereafter the behavior and main properties of Sb 
by assuming general arrays of water Cherenkov detectors. The results, how- 
ever, should be still qualitatively valid for any other detector which produces 
comparable signals for the muonic and electromagnetic components of the 
showers. Therefore, the signal and Sb are measured in VEM (Vertical Equiva- 
lent Muons, as in Haverah Park and Auger experiments). In this case, we will 
demonstrate later that the primary identity discrimination power goes through 
a maximum around b = 3. The analysis has been performed considering differ- 
ent array geometries and varying the detectors separation, showing the general 
applicability of the parameter. Sb could also be applied to scintillator arrays 
as the future Telescope Array experiment [5]. 

The paper is organized as follows: Section 2 shows an analytical discussion of 
the properties and stability of the new parameter. Section 3 shows in some 
respects an equivalent numerical study performed with simulations taking into 
account the effects of reconstruction and experimental uncertainties. In Section 
4 we perform a realistic comparative study among the stability and reliability 
of the inferred composition by using S 3 , X max and the rise time at 1000 m 
from core. Conclusions are presented in Section 5. 



2 Analytical study 

The parameter S& for a given event is constructed from the total signal in each 
triggered Cherenkov detector. Therefore, it depends on the normalization and 
shape of the lateral distribution function of the total signal. Close to the 
impact point of the shower, the signal is dominated by the electromagnetic 
particles (photons, electrons and positrons) whereas at larger distances it is 
dominated by muons. Fig. 1 shows the muon, electromagnetic and total signal 
in the Cherenkov detectors as a function of the distance to the shower axis for 
protons and Iron nuclei. The zenith angle of the simulated events considered is 
such that 1 < sec6 < 1.2 and the primary energy 19 < log(E/eV) < 19.1 (see 
section 3 for details about the simulations). The hadronic model considered is 
QGSJET-II. Fig. 1 also shows the fits of the LDF of each component with a 
NKG-like function [21,22,23], 
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where r s = 700 m and tq = 1000 m, and So and (3 are free fit parameters. 
If we consider proton and Iron primaries, the discrimination power of a mass 




r [in] r [m] 

Fig. 1. Lateral distribution functions of the muon, electromagnetic and total signal 
in the Cherenkov detectors for simulated protons (left) and iron nuclei (right) of 
1 < sec (9 < 1.2 and 19 < log(£'/eV) < 19.1. The hadronic interaction model used 
to generate the showers is QGSJET-II. The solid lines correspond to the fits with a 
NKG-like function (see Eq. (2)). 

sensitive parameter q, like Sb, can be estimated by using the so-called merit 
factor defined as 

n = , E[qfe] ~ SM , (3) 

where E[qj^\ and yar[g^] are the mean value and the variance, respectively, 
of the distribution function of the parameter q& with A = pr, fe. Note that 
an alternative definition for the merit factor makes use of the median instead 
of the mean value and, instead of the variance, cr| 8 [g] = [(g 84 — g 16 )/2] 2 , where 
gs4 and qie are the quantiles corresponding to 84% and 16% of probability, 
respectively. Now we use the definition as it is in Eq. (3) to make possible the 
analytical analysis. 



2.1 Optimization in water Cherenkov detectors arrays 



Assuming that the fluctuations of the total signal in a water Cherenkov de- 
tector are Gaussian, the distribution function for a given configuration of 
triggered stations is given by, 



P(st, ...,s N ;ri,.. .,r N ) 



f(n, ...,r N ) 
(2n)^Ul 1 cr[S(r l )} 



exp 



f {Sj - S(r z )f 
h 2 ^[S( ri )] 



(4) 

where rj is the distance to the shower axis of the zth station (the first station, 
7*1, is the closest one), Sfc) is the average LDF evaluated at r i; cxfS^rj)] = 
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(1 + e) [^(r^/VEM] 1 / 2 VEM [14,22] and f(r u ...,r N )is the distribution func- 
tion of the distance of the different stations to the shower axis. We have 
checked that for values of e of at least 5% there is no change in the results. 
Note that just two of the random variables {ri, . . . ,tjv} are independent, for 
instance, choosing n and r 2 (the first and second closest stations to shower 
axis) as the independent ones, we can write f(r±, . . . , rjv) = jfi,2(ri, r 2 )S(r 3 — 
r 3( r i ; r 2)) ■ ■ ■ 5(tn — rjv(ri, V2)), where 5(x) is the Dirac delta function. 

From Eqs. (1) and (4) we obtain the expressions for the mean value and the 
variance of S&, 
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where the variables have already been integrated and, 
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Here, fiirj) is the distribution function of the distance to the shower axis for 
the ith station and fij(ri,rj) is the distribution function of the distance to 
the shower axis of the ith and jth stations, 



fi,j{n, rj) = d^... dri_ x dr i+x . . . drj-idr j+1 . . . dr N f(r u ...,r N ). (9) 



In order to simplify the expressions for the mean and variance of S^, we perform 
the following approximations, 



E[g{rij\ = g(E[ri\), 



cov[g(ri),g{rj)\ = — 



dr 



cov[ri,rj}, 



(10) 
(11) 
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where g(r) = S(r)(r/r ) b . Thus, we finally get, 
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Fig. 2. Distance of the stations to shower axis for almost vertical showers in a 
triangular grid of 1.5 km of spacing. 
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We already have analytical expressions for the average LDFs of proton and 
Iron primaries obtained by fitting the simulated data (Fig. 1). The other in- 
gredients needed to calculate the mean value and the variance of S b are the 
mean values of the distance to the shower axis for the different stations and 
the covariance between all pairs of those random variables. We obtain these 
quantities from a simple Monte Carlo simulation: we uniformly distribute im- 
pact points in a triangular grid of 1.5 km of spacing, like the Auger array, and 
then, for each event of zenith angle such that sec^ = 1.1 and azimuthal angle 
uniformly distributed in [0,27r], we calculate the distance of each station to 
the shower axis. The result is shown in Fig. 2. From these distributions E[rj\ 
and cov[ri,rj] are easily determined. 

Finally, we could calculate the mean and the variance of 5&, and therefore, 
the merit factor. Fig. 3 shows the discrimination power r] as a function of b 
obtained under the mentioned assumptions and simplifications. We see that t] 
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2.2 Modifying the slope of the LDF 



We also study the discrimination power of Sb when the slope parameter /3 (see 
Eq. (2)) is modified but keeping constant the integrated signal for distances 
larger than the Moliere radius, ru = 80 m. Thus, the inferred energy of the 
event by surface experiments would not be significantly affected. The modified 
LDF that fulfills this condition can be written as, 

Qf m N(r M ,r ,r s ,(3 ) 

S(ri/ " = JV(r„,ro,r.,/3) (M) 

where, 

r 2+2/3 

N(r Ml r ,r s ,P) = Beta(-r,/r M , -2(1 + /3), 1 + /3), (15) 

ro{r s + r ) p 

Beta(*,a,6) = f*dt ^(1 - t) 6 " 1 and S 0o {r) is the LDF of Eq. (2) with the 
parameters So and (3q originally obtained from the fits in Fig. 1. 

The slope of the proton LDF is smaller than the corresponding to iron (the 
absolute value is greater). Then, we modify the slope of both LDFs such that, 
/V(0 = P%- ~ (£ - l)AA>/2 and /3 /e (0 = f3% + (£ - l)A/3 /2, where (3° pr and 
/3° e are the proton and Iron slopes, respectively, obtained from the fits of Fig. 
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Fig. 4. Contour plot of r/(£, b)/rj(l, 3). 

1, A/3 = /3/ e — fip r and £ is such that A/3(£) = £A/3 , i.e., £ = 1 corresponds 
to the non modified case. Note that for £ = 0, /3 pr = /3/ e = (/3 pr + /3/ e )/2. 
The mean and the variance of Sb are calculated using the same procedure 
as before, but the previous S(E[ri\) values are now modified by the factor 
N(r M , r , r„ /3 )/N(r M , r , r s , 0). 

Fig. 4 shows a contour plot of t/(£, b)/rj(l, 3) from where it can be seen that 
as £ increases 77 also increases. We also see that the maximum of 77 remains 
close to b = 3 almost independent of £. 



3 Modifying the muon content of the simulated showers 

There is experimental evidence of a deficit in the muon content of the simulated 
showers [24,19]. It is believed that such deficit is originated in the high energy 
hadronic interaction models which are extrapolations, over several orders of 
magnitude, of lower energy accelerator data. As mentioned, the total signal 
can be decomposed in the muon and electromagnetic signal. Therefore, to 
study how Sb chang function of the muon content of the showers we 

modify the total LDFs in the following way: S(r) = S em (r) + fS^(r) where / 
parametrizes the artificial variation in the muon component. 



9 



550 




Iron 



500 



450 — 



W 400- 

> - 



Proton 



W3 350 



300- 



250- 



IQ 

0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

/ 



Fig. 5. Mean values of S3 for protons and iron nuclei as a function of / where / = 1 
corresponds to the muon content predicted by QGSJET-II. 

The mean and the variance of Sb are calculated following the same approxi- 
mations explained before. For example, the mean value is given by, 



and similarly for the variance. The signal of the electromagnetic and the 
muonic components were also fitted separately, see Fig. 1. 

Fig. 5 shows the mean values of S3 for protons and Iron nuclei as a function 
of /. As expected, they increase with /. We also see that the iron curve 
increases faster than the proton one, which means that, for larger values of /, 
the discrimination power of £3 also increases. This happens because the muon 
content of the showers is very sensitive to the primary mass. Then, for larger 
values of / the muon component becomes more important increasing the mass 
sensitivity of S3. 

Fig. 6 shows a contour plot of rj(f,b)/r](l,3) from which we see how the dis- 
crimination power of Sb increases with the muon content of the showers and 
that the maximum is reached at b = 3 almost independently of /. 
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3 Numerical analysis 



3. 1 Sb under different array geometries 



The detection of the extensive air showers by a surface array of water Cherenkov 
tanks has been simulated by using our own simulation program previously re- 
ported in [25]. This program allows us to change easily the geometry of the 
array and the distance between detectors. Thus, triangular and square grids 
have been considered varying the array spacing (A) from 500 to 2000 meters. 
The injected LDF used to assign the signal to each detector is given by Eq. 
(2), where Sq and for both primaries are taken from the fits shown in Fig. 1. 
Thus, we assume here almost vertical showers with energy around 10 190 eV. 

We use again the merit factor, rj (see Eq. (1)), to quantify the discrimination 
power of Sb as a function of the exponent b but, contrary to the previous 
Section, we prefer now to use the median and <7g 8 instead of the mean value 
and the variance, respectively. The result for a triangular grid is shown in 
Fig. 7(a). It can be seen that, despite the fact that the value of r] decreases 
markedly as the density of the array diminishes, the maximum discrimination 
power is reached for b = 3 almost independently of the array spacing. The 
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Sb exponent Array Spacing [m] 

Fig. 7. Left: Merit factor as a function of Sb exponent for several distance between 
detectors in a triangular grid. Right: S& merit factor for the optimum b as a function 
of the array spacing for a triangular and a square grid. 

result is in agreement with the analytical approach for the 1500 m array (see 
Fig. 3), although the merit factor is now smaller because shower fluctuations 
are also taken into account in the present calculation. Similar results have 
been obtained in the case of a square grid, showing that S& encodes mainly 
information about the shower front structure and not just about the array 
geometry. 

The maximum merit factor as a function of the array spacing is shown in Fig. 
7(b) for both the triangular and square array. The merit factor increases as 
the array spacing decreases, as it was expected since the LDF is sampled in 
more points as the array becomes denser, improving Sb separation power. In 
case of A < 1000 m the merit factor is larger for the triangular grid since the 
number of triggered stations is also larger for this geometry. 

3.2 Simulations 

In what follows, we perform a more realistic simulations in order to treat 
more accurately the tank response and to take into account experimental 
uncertainties such as the shower reconstruction and the hadronic interaction 
model. In addition, we extend the energy and zenith angle range. 

The simulation of atmospheric showers is performed by using the AIRES 
Monte Carlo program (version 2.8.4a) [27] with QGSJET-II [28] and Sibyll 
2.1 [29] as the hadronic interaction models (HIM). Since the number of sec- 
ondary particles produced in a shower is extremely large (i.e. ~ 10 11 particles 
in a proton shower of 10 20 eV), it is very costly, in processing time and disk 
space, to follow all of them. Therefore, we use a statistical method called thin- 
ning, first introduced by M. Hillas [30], as it is implemented in AIRES. A 
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relative thinning of 10~ 6 and weight factor of 0.2 are used for the generation 
of the showers. Iron and proton primaries are simulated in the energy range 
from 10 19 to 10 19 ' 6 eV and the arrival direction follows an isotropic distribu- 
tion with zenith angle in the range 0° < 9 < 60°. The number of simulated 
showers per HIM and primary corresponds to an exposure of 5165 km 2 sryr 
(~ 0.8 years of Auger full operation [32]). 

The simulation of the response of the surface detectors as well as the shower re- 
construction are performed using the official Offline reconstruction framework 
of the Pierre Auger Observatory [33]. The simulation includes a triangular 
grid of Cherenkov detectors of 1.5 km of spacing. The unthinning method of 
P. Billoir [34] is used to obtain a list of particles that hit a given detector of the 
array. The GEANT4 package [35] is used to simulate the behavior of particles 
inside the tanks. The surface detector simulation has been tested and proved 
to be in good agreement with experimental data [36] . In order to increase the 
statistics, each shower is recycled 5 times by randomly distributing their cores 
inside the array (full discussion of the statistical effects of recycling air showers 
could be found in [37]). 

The reconstructed energy has been simulated by fluctuating the real one using 
a Gaussian uncertainty of 20%, a typical value for surface experiments [31,32]. 

Simulations are divided in logarithmic energy bins from log(E/eV) = 19 to 
\og(E/eV) = 19.6 in steps of 0.1. We also consider three different bins of 
zenith angles centered at 30°, 45°, and 55° and of 10° wide. 

3. 3 Optimization 

Fig. 8 shows the merit factor of as a function of b for both HIM considered. 
The maximum is reached at b = 3 in a very good agreement with the result 
obtained in the analytical analysis (see Fig. 3). Furthermore, the shape of 
the curve is quite similar. Nevertheless, the merit factor is lower and the 
peak wider, as it is expected because every event, independently of its energy 
or zenith angle, has been included and shower fluctuations and the effects 
introduced by the detectors and reconstruction methods are accounted for. 
The result is also in agreement with previous numerical simulation (Fig. 7(a)). 



3.4 Influence of distant detectors 

Several tests have been performed to study whether distant stations from 
shower axis, for which the fluctuations in the total signal are quite significant, 
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Fig. 8. Merit factor of Sf, as a function of b obtained from simulated events. There 
is a general agreement with Fig. 3 and Fig. 7(a) 

could affect the separation power of S3. Note that in the S3 sum, saturated 
stations are not included. 



Let us call ru m to the maximum distance of the stations considered in the 
S3 sum. In previous analysis, no cut in distance were applied so that all the 
triggered detectors were included with the same weight. Now, different cuts 
are tested: 

• r Hm = r o P t- We call r opt to the distance obtained by searching the value 
that maximizes the merit factor in 100 m steps. It depends on the hadronic 
interaction model and zenith angle. 

• fum = 2700m. All the stations at a distance from shower core larger than 
2700 m are excluded. This value is selected, somehow arbitrarily, in order 
to reject the stations whose signal could be dominated by fluctuations. 

• fu m = r max (E,9). Considering the distribution of the distance from the 
furthest triggered detector to the shower axis (see Fig. 9), r max is defined as 
the median of such distribution. In the calculation for each energy and zenith 
angle intervals, proton and Iron primaries and both HIM are included. 

• Each term of the S3 sum is weighted using the so called Lateral Trigger 
Probability (LTP), i.e. the probability of a shower to trigger a detector as a 
function of the reconstructed energy, shower zenith angle and the distance 
from the detector to the shower axis. This is an elegant way of switching off 
smoothly the stations at large distances from the core. The LTP is obtained 
as described in [38]. 

• fu m — » 00. All the triggered stations are included. No cut is applied. 

Fig. 10 shows r) as a function of the primary energy for the different cuts men- 
tioned before. It can be seen that the difference among them is not significant, 
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Fig. 9. Distance from the furthest triggered station to shower axis for three dif- 
ferent zenith angles as a function of energy. Both hadronic interaction models and 
primaries are included. The error bars correspond to the 68% and 95% confidence 
levels. 
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Fig. 10. Merit factor as function of the reconstructed energy. Several cuts are tested 
(see the text for details). Here the result for Sibyll 2.1 and 9 = 55 ± 5° is shown. 
Similar result are obtained for QGSJET-II and other zenith angle bins. 



showing the robustness of the parameter. Therefore, the discrimination power 
of S3 is not affected by distant stations where the fluctuations could dominate 
the signal. Since the cuts do not improve the parameter ability, we recommend 
no to use a cut in distance. Furthermore, these cuts could introduce bias in 
composition determination since the LDF of a proton primary is steeper than 
that of an Iron nuclei of the same energy and zenith angle (see Fig 1). 
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Fig. 11. S3 vs. sec{9) for proton primaries (Sibyll 2.1). A linear fit to the points 
shows a negligible dependence on the zenith angle of the incoming shower. Error 
bars are RMS/VW. 
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Fig. 12. log(S 3 /VEM) vs. log( J B/eF) for = 30° (left) and 9 = 45° (right). The 
error bars are the RMS/y/N- Linear fits to the points are also shown. 



3.5 S3 dependence on the primary energy and zenith angle 



Fig. 11 shows the mean value of S3 as a function of 360(6*) for protons primaries 
and Sibyll 2.1. There is no significant dependence with the zenith angle. Sim- 
ilarly occurs in case of Iron primaries and QGSJET-II. 



The evolution of S3 with energy for proton and iron primaries is shown in Fig. 
12 for two zenith angle bins, centered at 9 = 30° and 9 = 45° and 10° wide. 
An almost linear dependence has been found. 
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4 A realistic application 



The reliability of the composition determination by using 63 and other mass 
sensitive parameters is checked in this Section. We select the two parameters 
most commonly used in the literature, one from the surface detectors, the 
rise time at 1000 m from shower core, and the other from fluorescence tele- 
scopes, X max , the atmospheric depth at which the maximum development of 
the cascade is reached. A brief discussion follows on specific details about the 
determination of these two parameters: 

• Rise time at = 1000 m from core, ^1/2(^0) [ ns ] : the rise time of a single 
station is defined as the time it takes to increase from 10% to 50% of the 
total signal. The spread of the arrival times of the shower particles at a 
fixed core distance increases for smaller production heights, so the rise time 
is expected to be smaller for heavy primaries that develop higher in the 
atmosphere. £1/2(7*0) is obtained by fitting the rise time of each triggered 
station using the function £1/2(7") = (40 + ar + br 2 ) ns and evaluating the 
fitting function at r$ = 1000 m. Parameters a and b are free in the fit. Only 
stations whose distance to shower axis is in the range from 600 to 1500 m 
and whose signal is greater than 10 VEM are included in this fit in order 
to avoid signals dominated by large fluctuations and saturated ones. Since 
three stations are required at least in the fit, the number of events for which 
is possible to evaluate £1/2(7*0) is significantly reduced, specially at larger 
zenith angles. 

• X max [g/cm 2 ]: up to now, X max is regarded as the most useful parameter 
for composition analysis. X max is very sensitive to the identity of the pri- 
mary and less affected by uncertainties in energy determination compared 
to surface parameters, since it depends logarithmically on primary energy. 
In order to calculate a realistic X max , we use the value calculated by AIRES 
and fluctuate it by using a Gaussian distribution whose standard deviation 
is a typical value for the resolution achieved in fluorescence experiments 
[11,12], accounting for the response of the detector and the effects of the 
reconstruction method. 

The parameters S3 and X max are almost independent of zenith angle, so it 
is possible to combine events of different zenith angles in a given sample. 
Obviously, this is not the situation for t 1 / 2 (r Q ). To compensate this dependence 
a quadratic fit is performed £1/2(7*0) vs - sec(8) for each primary and hadronic 
model, and using average values of the fitted parameters, we correct, in a 
simple way, for the zenith angle dependence of £1/2(7*0): 

tT/ 2 r (r , sec 0) = tT/Tiro, sec 6) + (t{%{r , 1.05) - t{%(r , sec 9)) . (17) 



17 



The correction does not increase the fluctuations and ^i°2 r ( r o) shows a strong 
reduction on the zenith angle dependence. 

For the subsequent analysis, we consider the lowest energy bin (form 10 19 
to 10 19 ' 1 eV) where we have the largest number of events. The samples cor- 
responding to a given primary, HIM and one of the mass sensitive param- 
eters considered are binned. Let us call h p (i) and hf e (i) to the number of 
events in the ith bin, normalized to the total number of events in the sample, 
corresponding to a given parameter and for protons and iron nuclei, respec- 
tively. The histograms h p and hf e are assumed to be the distribution of the 
universe. The proton abundance or composition of a sample is defined as 
C p = N p /(N p + Nj e ), where N p and Nf e are the number of protons and iron 
nuclei in a given sample. We consider samples of Cp Ue from to 1 in steps of 
0.1. For each value of C* rue , we generate 500 samples of N s = 300 events each 
by taking at random values from the histograms h p and hf e . For each of these 
new samples we generate a histogram H s , with the same binning used in h p 
and hf e , which is not normalized and then satisfies J2iH s (i) = N s . Thus, as- 
suming Poisson statistics, the probability of having a configuration with H s (i) 
events in the ith bin is given by, 

W)W = II «*(-£*(«')) x ^fTTT ( l8 ) 

where H t {i) = N s (C p h p (i) + (1 — C p )hf e (i)) and C p is the unknown pa- 
rameter. Therefore, the inferred proton abundance is obtained by minimizing 
-In P({H s (i)}i). 

Fig. 13 shows the inferred composition as a function of the true composition 
corresponding to QGSJET-II. For Sibyll 2.1 similar results are obtained but 
with smaller error bars (which represent the 68% and 95% C.L.), in agreement 
with the fact that the merit factors are in general greater for this HIM. 

The statistics available for each mass sensitive parameter corresponding to a 
given exposure time, is a key value to compare their discrimination capabilities. 
Due to the limited duty cycle of the fluorescence telescopes, only 10% of the 
events are detected, so the statistics for X max is significantly lower than that 
for surface parameters (ground detectors have almost full operation time). 
This fact has been considered in Fig. 13. Just to illustrate the significance 
of taking into account the limited statistics for X max when doing composition 
studies, it is also shown the result for X max if the same statistics as the surface 
parameters were available. Then, the error bars are reduced becoming the 
smallest ones, but ten times more exposure would be required. 

A second study has been performed. Now, a fix true proton fraction Cp Ue = 0.5 
is assumed and the inferred proton fraction is calculated in the energy range 
from 10 19 to 10 19 ' 6 eV. In this case, in order to improve the small statistics 
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Fig. 13. True vs. inferred proton fraction at ~ 10 19 eV using QGSJET-II for ii/2( r o) 
(top-left), S3 (top-right) and X max (bottom- left) for a fixed exposure time (the 10% 
duty cycle of the fluorescence telescopes is taking into account). The bottom-right 
panel shows the inferred proton fraction obtained using X max but for samples 10 
times larger, corresponding to the same SD exposure. 

in the higher energy bins, the histograms h p (i) and hf e (i) corresponding to a 
given parameter, energy bin and HIM are fitted. The fitting function used is 
the so-called Asymmetric Generalized Gaussian (AGG) [39] which is defined 

as, 



Wv) = r wc ' Hz + aO'] «7y<A* (19) 




where, 

1 /r(3/c) 



1/2 



la 



Oi + O r \r(l/c) 




of and of are the variances of the left and right sides of the probability density 
function, respectively, and T(x) is the Gamma function. If of = of, the AGG 
is symmetric. Furthermore, if of = o^ and c = 2, AGG reduces to the reg- 
ular Gaussian distribution function and, for c = 1, it represents the Laplace 
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Fig. 14. Examples of the fits with AGG function for the three parameters considered 
for different energy bins and hadronic interaction models. 

distribution. 

Fig. 14 shows one example of these fits for each parameter considered. As can 
be seen, the AGG function allows us to fit asymmetric distributions with longer 
tails and sharper (or natter) maximums compared to a Gaussian distribution. 
By using the distribution functions obtained in these fits, we could get the 
needed events to perform the calculation with enough statistics. Samples of 
events corresponding to a given primary type, parameter, energy bin and HIM 
are generated by sampling the corresponding fitting function. The histograms 
hp and hf e used as the universe have 1000 events each. The number of events 
of the test samples (used to calculate the error of the inferred composition) 
varies as a function of primary energy because of the steepness of the spectrum. 
The number of events expected by Auger in 1 and 5 years of full operation 
are considered"*"! and to reproduce real conditions, the number of events with 
available X max is 10% of the total in the sample. The procedure to infer the 
composition is the same as explained before. 

As in previous case, there is no significant bias in the inferred proton abun- 
dance. Fig. 15 shows the uncertainty (at 68% of C.L.) on the determina- 
tion of the proton abundance as a function of primary energy for samples 



For example, in one year of full operation and considering the spectrum reported 
in [32], around 400 events are expected at 10 19 eV and around 30 at 10 19,6 eV 
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Fig. 15. Error in the inferred proton abundance determined by using S3, X max and 
ti/2( r o) f° r 1 and 5 years of Auger full operation time. Top: QGSJet-II. Down: 
Sibyll 2.1. The statistics used for X max is only 10% compared to that for surface 
parameters due to the limited duty cycle of the fluorescence telescopes. 

of Cp Tue = 0.5. It can be seen that the best results are obtained for S3. As 
mentioned before, the errors corresponding to Sibyll 2.1 are smaller than for 
QGSJET-II because the shower to shower fluctuations are in general smaller 
for Sibyll 2.1. 

As mentioned in Section 3.2, the reconstructed energy considered in this work 
is obtained fluctuating the real one with a Gaussian uncertainty. This proce- 
dure does not take into account the possible correlation between the recon- 
structed energy and primary composition, which has the potential to degrade 
the discrimination power of any mass composition parameter. In order to as- 
sess the magnitude of such effect, we consider the energy estimator S38 as in 
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Fig. 16. AS"38 = 5 < 38(/e) — S3g(p) (left) and AS^/ Sssip) (right) as a function of the 
real energy. The error bars are the RMS/ yN. 



the case of Auger [7]. 



It is expected that, for the same energy, different primaries will produce show- 
ers with different S(r) parameter, giving rise to systematic effects on the re- 
constructed energy depending on the unknown composition of the impinging 
cosmic ray flux. In fact, from the simulations, we calculate the difference in 
5*38 between iron and proton showers with the same real energy and, as Fig. 
16(a) shows, ASss = S-^fe) — 638 (p) increases almost linearly with the real 
primary energy. The relative value AS^/ S^ip) is also shown in Fig. 16(b) 
(which is in agreement with [40]). This difference translates into a spurious 
difference, AE, between the estimated energy for proton and Iron primaries 
of the same arrival energy. Therefore, in order to assess the influence of the 
correlation between the reconstructed energy and composition on the stability 
of our S3 parameter, we modify the reconstructed energy of each iron event by 
adding to its originally estimated energy the systematic error, AE, calculated 
from its own S38. On the other hand, the procedure for the reconstruction of 
the energy of proton primaries is not modified. The energies thus obtained are 
used to repeat the process described above for composition discrimination in 
the case of surface parameters. 



Fig. 17 shows the results obtained for 1 and 5 years of integration time. The 
error in the measured composition using S3 increases but still remains smaller, 
or of the order of those obtained using rise time or X max as discriminator 
parameters. The results shown in Fig. 17 were calculated using Sibyll 2.1 as 
hadronic interaction model, but are qualitatively similar for QGSJET-II. 
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Fig. 17. As Fig. 15(b) but considering the correlation between the reconstructed 
energy and the primary composition. See text for details. 

5 Conclusions 



A new family of parameters, that we call Sb, has been proposed for compo- 
sition analysis in cosmic ray surface arrays. The parameters are determined 
on a shower-to-shower basis from the total signal deposited in each triggered 
detector and their distance to the shower axis. Despite the fact that surface 
composition parameters are usually more affected by systematic errors than 
the corresponding parameters obtained from the fluorescence technique, the 
former are of great interest because of their larger duty cycle (around 10 times 
larger) . 

Sb is of general applicability: its discrimination power has been studied for 
different surface arrays, varying the distance between detectors and the array 
geometry (square and triangular grids). The separation power of Sb increases 
rapidly as the array spacing decreases and remains strong for any array ge- 
ometry. 

We have performed an extensive analytical study of the most relevant proper- 
ties of Sb- We demonstrate that there is a well defined value of the parameter 
b which maximizes the discrimination capability of Sb- In the case of a binary 
mixture or proton and Iron, that optimal value is b ~ 3. 

Furthermore, motivated by experimental results that suggest an excess of 
muon shower size with respect to numerical simulation expectations, we have 
analyzed the stability of b under such conditions, showing that the optimal fit 
for b is still ~ 3 and that, furthermore, the discrimination power of S3, for an 
Fe-p mixture, is actually enhanced. 
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Numerical simulations have also been carried out to account for experimental 
conditions and the uncertainties involved in the reconstruction procedure. The 
previous analytical results have also been validated from the simulations which 
include the detector response and reconstruction efficiency. The numerical 
results support that b = 3 is the value that maximizes the discrimination 
power of Sb, in agreement with the analytical method, and independently of 
the array geometry or the distance between detectors. 

Additionally, we demonstrate that 5*3 is almost independent on zenith angle, 
which represents an advantage over other mass sensitive parameters, and it is 
almost linearly dependent on the primary energy. 

A realistic comparison between 5*3 and the most useful parameters from the 
surface and fluorescence technique, i.e., the rise time, tyz, and X max respec- 
tively, has been performed. Our simulations show that the accuracy on the 
determination of the composition obtained by using S3 is, at constant integra- 
tion time, greater than the one obtained using the other parameters consid- 
ered. Even if a realistic correlation between the reconstructed energy and the 
primary composition is included in the analysis, the discrimination power of 
5*3, albeit degraded, still remains greater or of the order of the one obtained 
using the other parameters. 

The application of S3 or any other composition parameter in a real experi- 
ment is not straightforward. The most significant problems that could affect 
S3 have been considered in this work, such as the energy uncertainty, zenith 
angle dependence, surface array density and geometry, the effect of detec- 
tors very far from the core and the uncertainties in the hadronic interaction 
models (specially in the muon component). On the other hand, X max is less 
prone to systematic errors but several circumstances must be considered to 
achieve reliable results. For instance, the longitudinal profile reconstruction 
suffers significant uncertainties coming from the fluorescence yield or possible 
direct Cherenkov light towards the telescopes. Atmospheric monitoring is also 
needed to account for the effect of clouds and aerosols. In hybrid experiments 
it requires an optimum synchronization with a surface detector as well. In 
addition, only those events whose X max is in the field of view of the telescopes 
could be properly reconstructed and this selection cut could introduce artificial 
biases that must be carefully corrected [11]. 
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